Efficient Recursion Method for Inverting Overlap Matrix

A new O(N) algorithm based on a recursion method, in which the computational effort is proportional to the number of atoms N , is presented for calculating the inverse of an overlap matrix which is needed in electronic structure calculations with the the non-orthogonal localized basis set. This efficient inverting method can be incorporated in several O(N) methods for diagonalization of a generalized secular equation. By studying convergence properties of the 1-norm of an error matrix for diamond and fcc Al, this method is compared to three other O(N) methods (the divide method, Taylor expansion method, and Hotelling’s method) with regard to computational accuracy and efficiency within the density functional theory. The test calculations show that the new method is about one-hundred times faster than the divide method in computational time to achieve the same convergence for both diamond and fcc Al, while the Taylor expansion method and Hotelling’s method suffer from numerical instabilities in most cases. The development of O(N) methods [1–10] and the revival of localized orbitals as a basis set [11–19] have been made during the last decade in order to extend the applicability of the first-principles molecular dynamics (FPMD) simulations using the plane wave expansion 1 and the Car-Parrinello method within density functional theories (DFT) [20]. However, only few applications of these O(N) methods to large systems have been reported within the DFT calculations [18,21,22]. Although there are a lot of limitations of the method based on the localized description [22], one of the limitations is that several O(N) methods require evaluating the inverse of the overlap matrix S which comes from non-orthogonality among the localized orbitals. In the generalized Fermi operator expansion (FOE) method [4] to the non-orthogonal basis we need to calculate the inverse of overlap matrix to construct the modified Hamiltonian H ′ ≡ SH , while Stephan et al. have proposed solving a linear equation SH ′ = H with the cutoff radii of H instead of calculating the inverse of overlap matrix. In the density matrix (DM) method [8–10] which is a promising approach for materials with a wide gap, fortunately, the evaluation of the inverse is not required during the optimization of grand potentials, although we have to evaluate the inverse of the overlap matrix for a good initial guess of the density matrix [10]. The block bond-order potential (BOP) method [2], which has good convergence properties for both insulators and metals, also requires the evaluation of the modified Hamiltonian H ′ as in method the FOE method. If the overlap matrix is sparse, the computational cost scales as the second power of the number of atoms N in the inverse calculation. Therefore, an efficient O(N) method for inverting the overlap matrix should be developed. So far, several O(N) inverting methods have been proposed. Gibson et al. used a simple method in which a linear equation SH ′ = H constructed for a finite cluster is solved without explicit calculation of S [23]. Mauri et al. considered approximating the inverse of overlap matrix by the Taylor expansion [7]. The approach could be an O(N) inverting method when the matrix elements in the pth moment O of the overlap matrix O are cut at a finite distance. Palser and Manolopoulos proposed to evaluate the inverse by Hotelling’s method which is similar to the iterative purification algorithm of the DM method [10]. The iterative calculation can be performed in O(N) operations, provided that the cutoff of matrix elements at a finite distance is introduced in the product of two matrices. It is worth pointing out that 2 the ideas of these O(N) inverting methods are analogous to those of the O(N) methods for the diagonalization. The divide method by Gibson et al. [23], the Taylor expansion method [7], and Hotelling’s method [10] strategically and mathematically correspond to the divide and conquer method [5], the FOE method [3,4], and the DM method [8–10], respectively. Therefore, one may expect that these O(N) inverting methods may have the convergence properties for realistic materials similar to the O(N) methods for the diagonalization [24]. However, it remains to be seen whether the expectation is meaningful or not. In this paper we propose a new O(N) method for calculating the inverse of the overlap matrix which is based on a resolvent and the block Lanczos algorithm. The new method is compared with the other three methods in terms of the computational accuracy and efficiency. Thus, our aim of this paper is to clarify the applicability of these four O(N) inverting methods for realistic materials. The paper is organized as follows. In Sec. II we present the theory of a new O(N) inverting method based on a recursion method, and also summarize the three other O(N) inverting methods. In Sec. III we discuss the convergence properties of these four O(N) inverting methods for the diamond and fcc Al within the DFT calculations using the 1-norm of an error matrix which will be related to the error in the eigenvalues in this section. In Sec. IV we conclude with clear characterization of the four O(N) inverse methods.